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Abstract 

In this paper we report on results in the study of spatially homogeneous 
cosmological models with elastic matter. We show that the behavior of 
elastic solutions is fundamentally different from that of perfect fluid solu- 
tions already in the case of locally rotationally symmetric (LRS) Bianchi 
type I models; this is true even when the elastic material resembles a per- 
fect fluid very closely. In particular, the approach to the initial singularity 
is characterized by an intricate oscillatory behavior of the scale factors, 
while the future asymptotic behavior is described by isotropization rates 
that differ significantly from those of perfect fluids. 



1 Introduction 



In cosmology, the matter model that is used most frequently is that of a perfect 
fluid, usually with a linear equation of state. This choice is quite general in 
Friedmann-Robertson- Walker space-times, where by inheritance of symmetry 
the stress-energy tensor must have the algebraic form of the stress-energy of a 
perfect fluid. The case is different, however, for anisotropic space-times: Con- 
sidering perfect fluids is a restriction and might be misleading, since it is unclear 
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in general, how robust the results on the behavior of perfect fluid solutions are 
under a change of the matter model. For example, it was shown in [7] that the 
past asymptotic behavior of Bianchi type I cosmological models with collision- 
less matter is considerably different from that of perfect fluid models. A more 
systematic analysis of the problem in the case of Bianchi type I models with 
general anisotropic matter sources has been carried out in 0] . Most importantly, 
this difference concerns the behavior of solutions in relation to the Taub (flat 
LRS Kasner) solution. While there exist perfect fluid solutions that are asymp- 
totic to the Taub solution as the initial singularity is approached, this behavior 
does not appear for solutions with anisotropic matter. Results of this type have 
broad ramifications, since the dynamics of spatially homogeneous cosmologies 
are generally conjectured to be the building blocks of the asymptotic dynamics 
of generic cosmological models (i.e., models without symmetries), see [8] for a 
recent discussion. The Taub solution plays a crucial role in this context; essen- 
tial ideas like the basic concept of Kasner eras (cf. the BKL asymptotics [1]) 
and the advanced concept of decay rates of inhomogeneities for generic models 
are connected with the Taub solution. In fact, the Taub solution already plays 
a crucial role in the treatment of Mixmaster dynamics, see [15]. It is therefore 
essential to see in which respect the dynamics of solutions toward the singularity 
is dependent on the choice of matter model, and whether the role of the Taub 
solution changes in this context. A systematic analysis of this problem in the 
case of Bianchi type I cosmological models has been carried out in [3] . 

In this paper we widen the analysis of [3] of the dynamics of cosmologi- 
cal models with elastic matter; the elastic matter model we consider is taken 
from [5j [TD] . This matter model is particularly suitable for our aims, since it 
contains perfect fluids as a limiting case and thus allows to directly compare the 
behavior of anisotropic matter solutions (elastic solutions) with the behavior of 
perfect fluid solutions. In [3], it has been shown that the asymptotics toward 
the singularity of elastic solutions is fundamentally different from that of perfect 
fluids: The approach to the singularity is oscillatory, hence, in particular, there 
do not exist any solutions that approach the Taub solution in the asymptotic 
limit. In this paper we study the past and future asymptotic behavior of elastic 
models for the locally rotationally symmetric (LRS) Bianchi type I case in full 
detail. In particular, we present a detailed analysis of the oscillatory behavior 
of the scale factors that determine the metric, and we describe the dependence 
of the amplitude of the oscillations on the properties of the elastic material; 
see Section [5] (As we show in [4], oscillations in the past asymptotics of Bianchi 
type I solutions are due to an 'overcritical' violation of the energy conditions in 
a neighbourhood of the initial singularity; we refer to the concluding remarks.) 

In addition to the dynamics toward the initial singularity we consider the 
dynamics of elastic cosmologies in the low density regime, i.e., the future asymp- 
totic behavior of models. Also in this context we observe a fundamental differ- 
ence between elastic solutions and perfect fluid solutions: Although isotropiza- 
tion occurs for all solutions, the isotropization rates of elastic models differ 
considerably from those of perfect fluid solutions. This remains true even in the 
case when the anisotropy properties of the elastic material are small and the 
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elastic matter thus resembles a perfect fluid very closely; see Section S) 

The methods we use in this paper are methods from the theory of dynamical 
systems, see, e.g., [T3]. In Section[3]we begin by briefly discussing the dynamical 
systems formulation for Bianchi type I elastic space-times. This approach is 
described in [3] in more detail; here, however, we adapt the formulation to our 
present purposes. Throughout the paper we use unit such that 8nG = 1 and 
c = l. 



2 Set-up 

Locally rotationally symmetric (LRS) Bianchi type I space-times are represented 
by a line element of the form 

- dt 2 + A 2 dx 2 + B 2 (dy 2 + dz 2 ) ; (1) 

we denote the spatial part of the metric by g. 

In the vacuum case, the solutions of the Einstein equations are the Kasner 
solutions 

- dt 2 +t 2a dx 2 +t 2b (dy 2 + dz 2 ) , (2) 

where a + 2b = 1 = a 2 + 2b 2 . There exist two different solutions, the non-flat 
LRS Kasner solution (a,b) = (—1/3,2/3), and the Taub solution (a, b) = (1,0). 
Since the Taub solution is a representation of a subset of Minkowski space- 
time, the Taub metric admits a (smooth) extension beyond t — to Minkowski 
space-time. In particular, the hypersurface t = is null and not spacelike. 

For perfect fluids with a linear equation of state p/ p = w three classes of 
solutions exist. The Friedmann-Robertson- Walker (FRW) solution is isotropic, 

- dt 2 + f V(3H-H»]) ( dx 2 + dy 2 + ^2) _ (3) 

The non-isotropic solutions isotropize and approach ([3]) as t — > oo (this will be 
discussed in detail below) ; toward the singularity these solutions are asymptot- 
ically vacuum, i.e., they approach Q as t — * 0. Solutions with A/A < B/B are 
asymptotic to the non-flat LRS Kasner solution; solutions with A/ A > B/B 
approach the Taub solution, 



A oc t 



l-{l + w)t 1 - w , B oc 1 + (2 - w)t 1 ~ w (4) 



as t — > 0. (The proof of (JH) is straightforward, when one uses the formalism 
introduced below.) These solutions possess a weak null singularity, i.e., like 
the Taub solution they admit a (continuous) extension of the metric beyond 
t = [15]. (However, p diverges as t — > 0.) 

In this paper we consider an LRS Bianchi type I space-time whose matter 
content is described by an anisotropic stress-energy tensor of the same algebraic 
type, i.e., 

T*i = diag (p, p a ,pb,Pb) ■ (5) 
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It is common to define the normalized principal pressures, 

Pa Pb 

w A = — , w)b = — , 
P P 

the isotropic pressure p, and 

p lp A + 2p B 1, . 

W = - = = - (WA + 2wb) ■ 

pip 3 

For perfect fluids, wa — wb = to; if this is not the case the stress-energy tensor is 
often said to describe an anisotropic fluid. However, the required specification 
of the principal pressures typically lacks a physical foundation. A choice of 
principal pressures is usually made ad hoc to simplify the Einstein equations. 

2.1 Elastic matter 

In this paper we shall consider a stress-energy tensor of the form (JSJ) that rep- 
resents elastic matter. The general relativistic theory of elasticity has been 
originally formulated by Carter and Quintana in [5] and further elaborated by 
by Kijowski/Magli [9], Beig/Schmidt [2j and Karlovini/Samuelsson [10 . Rel- 
ativistic elasticity is used in both relativistic astrophysics and cosmology, see, 
e.g., [10] and p]. 

An elastic material is specified by a constitutive equation (Lagrangian) that 
depends on scalar functions constructed from the configuration map between the 
space-time and the natural (unstrained) state of the material; the stress-energy 
tensor is then obtained as the variation w.r.t. the space-time metric of the 
matter action. In particular, this yields expressions for the principal pressures 
without the need to resort to any ad hoc assumptions; we refer to [10] and to the 
appendix at the end of this paper. Note in this context that the anisotropics 
of the elastic stress-energy tensor are not due to intrinsic anisotropics in 
the elastic matter model, but to anisotropics of the space-time (provided the 
unstrained state of the elastic material is assumed to be isotropic) . More specif- 
ically, for the constitutive equation used in [2 [9] the energy density is 

P = M ^)-u + »>( 1 + ^<4^|!)!), (6 ) 

where po > is a constant, and the (normalized) principal pressures are given 

by 

vw A 4 - B A 

vw A i -B i 

W B - W+— A2B2 + i ™ (A2 _ B2)2 • (7b) 

These relations contain two constants, v and w, where 

| to | < 1 , vw > 
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is assumed. In the appendix of this paper we give a complete derivation of the 
energy density ([6]) and the principal pressures ([7]) for our particular choice of 
elastic matter and initial data. 

The constant v we call elastic constant. It measures the response of the 
elastic material to anisotropies; when v = 0, the elastic material reduces to 
a perfect fluid with equation of state p = wp. Therefore, the elastic material 
under consideration contains, as a special case, the perfect fluid model most 
widely used in cosmology, see [IB] , 

By I©-©, there are two regimes in which elastic matter can be viewed as 
a small perturbation of a perfect fluid: When vw <C 1 (small shear) or when 
A ~ B (almost isotropic geometry). This observation implies that there are 
indeed regimes where the energy distribution of a perfect fluid and of elastic 
matter are effectively indiscernible. Yet, as we will see, these models give rise 
to a rather different dynamical behavior of the space-time metric in the limits 
toward the initial singularity and for late times. 



2.2 Dynamical systems formulation 

The Einstein equations for a metric of the type (TT]) with stress-energy tensor ([5]) 
form a system of ODEs for the variables (A, A, B, B); in addition there is one 
constraint equation. It is preferable, however, to replace (A, A, B,B) by scale- 
invariant variables to obtain a dynamical systems formulation of the equations 
on a compact state-space. To this end we note that Jg = AB 2 , and we define 
the Hubble scalar as 



We introduce the dimcnsionlcss variables 



B 
B 



(8) 



3A 3B 



(9) 



which satisfy the identity a + 2b = 1. In the vacuum case, (a, b) coincide with 
the constant Kasner parameters and the metric is ©; note, however, that in 
the general case, (a, 6) neither satisfy the relation a 2 + 2b 2 = 1, nor does the 
metric assume the form ([2]). In the variables (a, b) the Hamiltonian constraint 
equation reads 



p 



3H 2 2 

using that > and a + 2b = 1 we infer 



(1 - a 2 - 2b 2 ) ■ 



(10) 



ae(-§,i), 6e(0,|) 

as well as fi < 1; in addition, $1 = 1 iff a = 1/3 = 6, which is the case iff 
A/A = B/B, i.e., for the FRW solution. 
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As the second set of dimensionless variables we use 

A 2 B 2 

1 P 



(11) 



A 2 + B 2 ' ' A 2 +B 2 ' 
where a + f3 = 1 . Finally we introduce a new time variable r by defining 

d T = H- X d t (12) 

and we denote by a prime the differentiation w.r.t. r. 

Using these variables the Einstein equations decouple into the dimensional 
equation for the Hubble variable 



H' = -3H 



i-§(i-«0 



and a reduced system of dimensionless equations 



a' = 9a(l — a) (a — — 



a' = —Q 



(1 -w)(a- ^ - (w A -w) 



where is determined by the constraint (fTQ|) . 

9 / 1 

n = -(l-a)(a+- 

and wa is given by (JT]) expressed in the new variables, 

1 - 2a 



wa = w + 



3 a(l-a) + if (l-2a) 2 



The phase space associated with the dynamical system (fl4l) is 



C 



(-i l) x (0,1) 9 (a, a); 



(13) 

(14a) 
(14b) 



(15) 



since (THJ admits a smooth extension to £, it is beneficial to include the bound- 
ary dC in our analysis. 



3 Dynamical systems analysis 

The boundary dC of the state space can be represented as a rectangle. The four 
vertices are fixed points of the dynamical system (|14[k we denote these points 

by 

Qi:(o,a) = (-|,l) Ti : (a,a) = (1,1) (16a) 

Qo : (a, a) = (-|, 0) T : (a, a) = (1, 0) . (16b) 
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The four sides of the rectangle (where we exclude the fixed points) are invariant 
subsets. When a = —1/3 wc find a' < 0; a — 1 entails a' > 0; for a — 1 we find 
by using (TT5)) that 



-£l(l-w) 



l + w/3 



1 — w 



<0; 



analogously, a' > when a = 0. It follows that the flow on dC forms a hetero- 
clinic cycle: 

Qi < Ti 



(17) 



Qo 



Tn 



The solutions associated with the fixed points To, Ti, and the orbit To — * 
Ti can be interpreted as the Taub solution. This is because a = 1 (b = 0) 
yields H cx I /{it) via (fT2|) and (fT3|) . and accordingly, (J9]) leads to A oc t and 
£? = const, i.e., to the Taub solution. Analogously, the fixed points Qo, Qi, 
and the orbit Qo <— Qi are representations of the non-flat LRS Kasner solution, 
since a = —1/3 (and thus b = 2/3). 

In the interior of C there is one fixed point: 

F:(a,a) = (§,§). 

The solution associated with F is isotropic, since a — 1/3 (and thus b = 1/3, 
fi = 1); a = j3 = 1/2 entails wa — wb = W. Accordingly, the fixed point F 
represents the FRW perfect fluid solution ([3]) , which is associated with a perfect 
fluid with equation of state p = wp. 

To analyze the global dynamics of the system (TH)) on L we introduce the 
function 



1\-! 

3 



1 



vw (1 - 2a) z 



Z= (l-a)- 1 !- 

V ; V 3J [ 6 a(l-a) 

which is positive on £; in fact, Z = 1 = inf £ Z at F, Z > 1 on £\{F}, and 
sup£ Z = +oo = Z\qc- A straightforward calculation shows that Z is decreasing 
along all orbits (different from F), i.e., Z' < on £\{F}. The monotonicity 
principle [6l [16] implies that (i) each orbit in C converges to F as r — > oo; 
(ii) each orbit (different from F) approaches dC as r — > — oo. Since dC is 
given by (fT7|) , the past asymptotic behavior of solutions is represented by this 
heteroclinic cycle. This behavior of solutions is depicted in Figure [TJ 

Equations (Tr2"|) and (TT3"|) allow us to translate between r-time and syn- 
chronous time t. Since — 3H < H' < —3(1 + w)H/2, Equation (|T2|l can be 
integrated to yield a positive function t{r) that satisfies t \ as r — > — oo. 

Thus the interpretation of the results on the global dynamics of solutions is 
the following: Each LRS Bianchi type I model with elastic matter isotropizes 
toward the future (i.e., as t — > oo) and (to first order) behaves like an (infinitely 
diluted) isotropic perfect fluid solution in the asymptotic regime. (At higher 
orders, however, when we consider isotropization rates, the behavior of clastic 
solutions differs significantly from that of perfect fluid solutions; see below.) 
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Figure 1: A generic orbit in the phase space C. It is assumed that vw > 
(3/32)(l — w) 2 . Under this assumption we observe oscillatory approach toward 
the FRW solution represented by F. The dominant energy condition is violated 
in the shadowed region. 



Toward the singularity, i.e., as t — > 0, we observe oscillatory behavior of 
elastic cosmologies. Elastic models do not converge to either the Taub solution 
or the non-flat LRS solution, but exhibit oscillations between the two, which is 
a direct consequence of the approach of solutions to the heteroclinic cycle (fl~7|) . 

In the following we investigate the future and past asymptotics of elastic 
models in detail. 



4 Future asymptotics and isotropization rates 

The analysis of the future asymptotic dynamics is based on an investigation 
of the dynamical system (|14|) in the neighborhood of the equilibrium point F. 
The linearization of (fl"4|) at F possesses the eigenvalues Ai, A2 and associated 
eigenvectors Vi, V2, 



Ai,: 



-(l-w)T\ (1-w) 



32 



Vl,2 = -7T > 77 



3 

16 



(1 — w) =F \j (1 — w) 2 — vw 



It is immediate that, if vw ^ ^(1 — w) 2 , the eigenvalues of the linearization 
of the system at F are real and negative; in this case F is a stable node. If 
vw > ^7j(l — w) 2 the eigenvalues are complex (with negative real part); in this 
case the fixed point F is a stable focus (stable spiral) and the solutions' approach 
to F as t — > 00 is oscillatory; see Figure [TJ The late time behavior of clastic 
cosmologies is thus characterized by 
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(i) monotonic isotropization if vw ^ (3/32)(l — w) 2 ; 

(ii) oscillatory isotropization if vw > (3/32)(l — w) 2 . 

To compare the behavior of elastic models with the behavior of the associated 
perfect fluid solutions we consider elastic matter that behaves "almost like a 
perfect fluid", i.e., we choose v to be small, 

vw<t:(l-w) 2 . (18) 

Evidently, this assumption implies monotonic isotropization of solutions. It 
follows from (TT5T) that 



3 . . vw 

Xi = --(l-w), A 2 = -4- 

2 1 — w 



Ul-w) 



vw 



in the lowest order approximation. Accordingly, the solutions of the dynamical 
system in the neighborhood of F are approximately given by 



f J + ci«i e-l Al l T + c 2 v 2 e~^ T . (19) 

Consider first the solution determined by c\ ^ 0, c 2 — 0. This solution coin- 
cides with the general perfect fluid solution, i.e., the solution to the (decoupled) 
equations 



which are obtained from (TH)) by setting = wb = w and thus form the 
system of equations for LRS Bianchi type I perfect fluid models (satisfying the 
equation of state p = wp) . Integrating (|13p we obtain 



H oc 



e"f ( 1+ ™) T (l + const e -3(i-«)T^ f ( 20 ) 



which is subsequently used in Equation (TT2]) to yield the functional relation 
between r-time and synchronous time t: 

toce i(l+w)r f 1 + congt 6 -3(1- w )t\ ( 21 ) 



(In the case w = 1/3 the lower order term is of a different form.) 

Finally, using (^/g)' = Sy^, which follows from ([5]) and (TT2"|) . we obtain 



V5 = V^t~ (l + const i~ 2 T+^ j , (22a) 
^(*) = (V^o) 173 ^ 1 ^ 7 f 1 - 2 const i-^) , (22b) 
B{t) = (Vffo) 173 ^ 1 ^ 7 (l + constt"^) . (22c) 
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However, it is straightforward to see from (fTO|) that this solution is not generic 
among LRS Bianchi type I solutions with elastic matter. Since |Ai| ^> |A 2 |, 
in (|19p . the term ciDie~' Al ' r can be neglected with respect to C2«2e - ' A2 ' r , and 
thus the generic solution behaves asymptotically like a solution given by c\ = 0, 
ri : 0. 

Proceeding in complete analogy to above we find 



v 7 ? = V9ot^ + const , (23a) 

A(t) = (y%) 1/3 ti^ (l - 2 const t~*^) , (23b) 
B(t) = (Vffo) 173 ^ 1 ^ 7 (l + const t~%T^} , (23c) 

for the generic elastic model. Comparing (|22f and ([23]) we observe a fundamental 
difference in the isotropization rates which are defined as 

isow = (vsr 1/3 |^-u/s) 1/3 

(I) There exists one single solution to the LRS Bianchi type I equations with 
clastic matter that behaves (to second order) like a perfect fluid as t — > oo. 
In particular, the isotropization rates are identical: 

iso(i) cx t~ isoa = t-&* . (24) 



(II) For the generic solution the isotropization rate is not the one of a perfect 
fluid solution, but 

iso(f) cx ^ is ° rf = f"*^ . (25) 

In particular, since iso c i <C Isoh, isotropization occurs at a much slower rate. 

We see that, therefore, the generic late-time behavior of (non-isotropic) so- 
lutions with elastic matter is considerably different from the behavior of (non- 
isotropic) perfect fluid solutions. Since this is despite the fact that the material 
properties of the elastic matter resemble those of a perfect fluid, this result is 
interesting. 

In astrophysical applications, such as the modeling of relativistic stars, re- 
placing perfect fluid matter by elastic matter yields a family of models that 
includes the perfect fluid solutions as special cases (which are obtained by sim- 
ply letting the elastic constants go to zero, see [10]). Here we see that, in con- 
trast, in (spatially homogeneous) cosmology, elastic models behave qualitatively 
different from fluid models. 



5 Past asymptotics 

In the asymptotic regime r — > — oo (i.e., t — > 0), every solution of (fl4|) is 
described by the heteroclinic cycle (Tl"T|) to an increasing degree of accuracy. 



10 



Accordingly, we observe alternations between Taub phases and non-flat LRS 
phases: In Taub phases, which are associated with the orbit in Figure Q] being 
close to either of the fixed points To, Ti or the orbit To — > Ti connecting these 
points, the solution is approximated by a Taub solution, i.e., 

A oc t , B oc const ; (26a) 

in non-flat LRS phases, which are associated with the orbit being close to either 
of the fixed points Qo, Qi or the orbit Qo <— Qi, the solution is approximated 
by 

A oc r 1/3 , B cx t 2/3 . (26b) 

Elastic cosmological models will oscillate between phases (|26a[) and (|26b[) with 
a rapidly increasing frequency as t — > 0. This is a simple consequence of our 
considerations. 

A priori, oscillations between phases of the type (]26a|) and (|26b[) are compat- 
ible with a large variety of scenarios: There could exist any constants < c\ < 
C2 < oo such that liminf t ^o^4 = c\ and limsup^o^ = ci- In the following, 
however, we show that 



(i) lim A = , if w > ; 
t— »n 

. . . (27) 

(ii) liminf A = 0, lim sup A = oo , if w < . 

t- * t->-o 

Accordingly, the amplitude of the oscillations converges to zero if w > 0, and 
diverges if w < 0. Note that convergence of the scale factor ^4 in the case w > 
occurs despite the fact that there exist phases (|26b[) where ^4 is increasing. 

In order to prove the asserted asymptotic properties of the scale factor A, 
we analyze in detail the asymptotic behavior of orbits in C Let e > 0; we define 
an e-ncighborhood of a point (a, d) £ C as the set of all points (a, a) such that 
|a— d| < e and |a— a| < e. If e is sufficiently small, then the flow of the dynamical 
system (|14l) in an e-neighborhood of a fixed point can be approximated by the 
linearized system. 

Consider an arbitrary orbit 7 in C. To facilitate matters we invert the di- 
rection of time by defining f = — r; hence the orbit 7(f) = (a, a)(f) approaches 
the heteroclinic cycle (fT7"l) as f — * 00 (instead of r — > —00). There exists an in- 
creasing sequence of times f n , neN, such that 7(17) = (a, a){f n ) — (1 — e, a n ), 
where a n < e for all n. Our proximate aim is to analyze the sequences 

(t„)„ s n and (a„) ne N ■ (28) 

By construction, at time t„, the orbit 7 enters the e-neighborhood of the fixed 
point To. Using the linearized flow it is straightforward to compute that 7 
leaves this neighborhood again at time f n + Af, where Af = (1/6) log(e/a n ); 
furthermore, 7(f n + Af) is given by (l — e~( 1 ~ w ' ) / 2 an w ^ 2 ,e). Subsequently, 
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the orbit 7 is approximated by a linear perturbation of the boundary orbit 
T -» Tx. 

The straightforward calculations show that the passage from the e-neighbor- 
hood of To to the e-neighborhood of Ti takes the time Af = — (1/3) log e and 
that the orbit 7 enters the e-neighborhood of Ti at a point with coordinates 
(l — ct e~ ( ~ 1 ~ w ^ 2 a.n^ ', 1 — e) where ct is a positive constant that depends 
on the choice of e only. Proceeding in this way, i.e., by tracking the orbit 7 
through the e-neighborhoods of the fixed points and along the boundary orbits, 
we obtain straightforwardly that 



X-n+l 



Ce 



^ ( 3-.) /(1+ „) 



where C is a positive constant that depends on e (but is independent of n) . As 
a alternative to (|29|) we can write 

(-^J = (-f) ' (30) 

where e is a positive constant (independent of n) . Iteration yields 

(f ) = (f) ■ (31) 

Note that by choosing < e we achieve oco/e < 1. Similarly, for f n+1 — f n we 
find 

a a. 

log + const , (32) 



Tn+l 



1 3(l + w) : 
and therefore 

f„ - to = . log ^ (1 + 0(nC-")) , (33) 

3(1 — w) e v 7 

where C > 1, or, approximately, 

a„ = const e -3(i-«0*"» . (34a) 

In a completely analogous manner one can find a different sequence of times, 
which we denote by f n , and a associated sequence 7(7^) = (1 — e, a' n ), where 
a' n = a(r^), such that 



1 - const e - 3{1 -' w) ^ . (34b) 



By construction, a n is related to the minimum of a in the interval [f n _i, t„] 
(which corresponds to one oscillation). To see that, note that the minimum is 
attained somewhere along the orbit Qo — * To; using again the linearized flow 
in the neighborhood of that boundary orbit, we infer that min f eff n _ 1 ,r„] a(f) = 
ka n for some constant k independent of n. Analogously, a' n is related to the 
maximum of a, i.e., max fg [ f ' ^] a(f) = fc'a^ for some constant A;'. 
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Since 



.4 3 



V~9 



a 



(35) 



1 - a 



by (IT]), and ^/g = y^e~ 3f by © and (|T2]). it follows that 
A n = A(f n ) oc exp [-(2 - w)f n ] , 
A ' n = Mt'u) oc exp [-wf' n ] . 



(36a) 
(36b) 



Consider the case w > 0. Since is a measure for the maximum that A attains 
during one oscillation period (and A n a measure for the minimum) , we find that 
A(f) — > as f — * oo; in fact, (|36|) can be condensed into the statement that 



as f — » oo (or r — > — oo, i — > 0). Note incidentally that e _T oc t 1 ^ 3 as i — > 0; 
this can be shown via (|13p and a line of reasoning analogous to the above. We 
have thus proved the claim that, in the case w > 0, the scale factor A converges 
to zero; this convergence occurs despite the existence of phases (|26bj) where A 
is increasing. In the case w < 0, the conclusion from ([36| is that the amplitudes 
of the oscillations grow unboundedly in such a way that liminf t ^o A = and 
hmsupj_, A — oo. This concludes the proof of the claim (|2"T|) . 

Interestingly enough, the behavior of the scale factors is largely determined 
by the constant w, while the dependence on the elastic constant v is minor: It 
is hidden in the constants appearing in (|34|) and in the derived formulas (I36|) 
and (|37|) . To explain the quasi-independence of the qualitative asymptotic be- 
havior of the scale factors of the elastic constant v, we note that, in the asymp- 
totic regime, the orbit spends a rapidly increasing amount of r-time in the 
neighborhood of the fixed points, while the time that elapses during the passage 
from one fixed point to the other is fixed. Since the fixed points thus dominate 
the asymptotic evolution of solutions and since the flow in the neighborhood 
of the fixed points is independent of v, it is only the isotropic constant w that 
enters in the description the qualitative asymptotic behavior of the scale factors. 

Summarizing we see that the past asymptotic behavior of elastic cosmologies, 
as described by (|36|) and (|37|) , is in stark contrast to the behavior of perfect fluid 
solutions, which converge to either the Taub solution (when a > b) or to the 
non-flat LRS solution (when a < b) as t — > 0; see (j4]). The structure of the 
singularity is therefore completely different for elastic models. In particular, 
there does not exist the option of a weak null singularity. 

6 Discussion and conclusions 

In this paper we analyze locally rotationally symmetric (LRS) models of Bianchi 
type I with elastic matter. Since the elastic matter model naturally contains 
perfect fluid matter as a limiting case — the latter being the matter model most 
commonly used in cosmology — , we are able to directly compare the behavior 
of elastic models with the behavior of perfect fluid cosmologies. 



A(t) < const e 



— wr 



(37) 
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Toward the future all elastic models approach an infinitely dilute isotropic 
state. The approach to this state is oscillatory when the elastic constant is 
sufficiently large; when the constitutive equation of state of the elastic matter 
does not deviate considerably from that of a perfect fluid (i.e., when the modulus 
of rigidity is small) isotropization is monotonic. Interestingly enough, even in the 
latter case, isotropization occurs at a rate that is fundamentally different from 
the isotropization rates observed for perfect fluids models. This is in contrast 
to standard astrophysical applications where elastic materials produce models 
that closely display the perfect fluid behavior [TO]. In particular perfect fluid 
solutions are recovered by letting the elastic constant v — ► 0. In the cosmological 
context this is no longer true; as shown in this paper, the isotropization rates of 
elastic solutions are qualitatively different from those of perfect fluid solutions 

The past asymptotics of elastic models is also interesting. While perfect fluid 
solutions converge to either the Taub solution or the non-flat LRS solution, clas- 
tic models oscillate between these two states. Oscillatory behavior toward the 
initial singularity is well-known in the context of certain Bianchi cosmologies, 
most notably in Bianchi types VIII and IX. This asymptotic behavior is usually 
referred to as the Mixmaster behavior. However, Mixmaster oscillations are ab- 
sent when LRS symmetry is imposed: LRS solutions simply approach the Taub 
solution or the non-flat LRS solution. This is in stark contrast to behavior of 
elastic models analyzed in the present paper. These models exhibit oscillations 
despite the fact that they are LRS models. 

Elsewhere we show in more detail that there is no evident connection be- 
tween elastic oscillations and Mixmaster oscillations [3]. (It might be justified 
to say that the former are caused by the matter and the latter by geometry) . In 
the paper [3j we consider the case of generic (in fact, diagonal) Bianchi type I 
models, and we find an intricate network of oscillations that determine the past 
asymptotic behavior of these elastic cosmologies. Again, this chaotic oscilla- 
tory approach to the singularity is fundamentally different from the Mixmaster 
behavior. 

Despite the fact that elastic oscillations are not directly connected with 
Mixmaster oscillations, there might exist interesting consequences when one 
considers more general elastic models than Bianchi type I. Consider, e.g., elas- 
tic models of Bianchi type VIII or IX, say. In this case we expect Mixmaster 
behavior, which is characterized by oscillations between epochs where the be- 
havior is close to the behavior of Bianchi type I models. Since already the 
Bianchi type I models are oscillatory, as shown here and in [3], we might be 
confronted with a hierarchy of oscillations, i.e., oscillations between oscillations, 
where Mixmaster oscillations connect epochs of elastic oscillations. Whether 
this is indeed true remains to be investigated. 

Although it is not the purpose of this paper to propose elastic matter as a 
physically realistic matter model for the universe, we would like to conclude this 
paper by commenting on the 'physics' of elastic matter and elastic cosmologies. 
We have seen that the asymptotic behavior toward the initial singularity of 
elastic cosmological models differs qualitatively (and not merely quantitatively) 
from that of perfect fluid models. This is hardly surprising, since the differences 
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in the material properties of elastic matter and perfect fluid matter are most 
significant under extreme conditions, such as those found in the proximity of 
a singularity. One could expect that a description of matter as a solid elastic 
body might not be completely unrealistic under these conditions — elastic mat- 
ter is also used in the modeling of compact stellar objects like neutron stars. 
However, we note that the dynamics toward the singularity of the elastic cos- 
mologies under study are connected with the violation of energy conditions. 
The dominant energy condition reads \wa\ < 1 and \wb\ < 1- This condition 
is violated in a neighborhood of the boundaries a = and a = 1 of the state 
space; to see this we observe that 

wa I n = w + 2 > 1 , wa I n = w — 2 < — 1 . 

la— la=l 

Hence, by continuity, there is a neighborhood of the sets a = and a = 1 in 
the phase space where the dominant energy condition is violated, see the shad- 
owed region in Figure [TJ when v — > 0, this region becomes smaller and smaller 
(and eventually reduces to the empty set for v = 0). Note that we only have 
directional dominant energy condition violation, i.e., only one of the normalized 
principal pressures (either wa or wb) is bigger than one. In particular, the 
isotropic pressure, w — 1/3(wa + 2ws) always satisfies \w\ < 1. The violation 
of energy conditions is thus much milder than in the context of phantom fields 
or similar models. In contrast to energy condition violation toward the singu- 
larity, at late times the dominant energy condition is always satisfied. In fact, 
for every solution, the energy conditions are satisfied for all times larger than 
some given time. (For a sufficiently small elastic constant v, the energy condi- 
tions are satisfied already after Planck's time.) Finally, let us draw the reader's 
attention to [4], where a close connection between energy condition violation 
and oscillatory singularities is established. 

In general, it is unrealistic to expect that the description of the matter as an 
elastic material — represented by the constitutive equation of state (JSj) — remains 
true under extreme stresses. Under extreme conditions the material will loose 
its elastic properties and its behavior might deviate considerably from the one 
described. For instance, the assumption of a constitutive equation of state of the 
quasi-Hookean form is typically justified only under the condition of small shear, 
see [TU], and hence violated if the shear scalar s 2 , cf. (j¥T|) . becomes too large; 
this occurs when the variable a is too close to or 1 (which characterises for the 
approach to the singularity). It is thus to be expected that the description of 
a material as elastic has a limited range of validity — it breaks down before the 
singularity is reached. Clearly, a modification of the material's properties under 
extreme conditions will lead to different asymptotic behavior of the associated 
cosmological models. We refrain from investigating these issues further here 
since any modification of the matter properties seems rather ad hoc rather than 
based on sound physical considerations. However, we refer to [4], where we 
consider more general anisotropic matter sources. 

The observed results on the long-term evolution of elastic cosmological mod- 
els are completely independent of the above considerations, especially since we 
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study elastic matter whose material properties are almost indistinguishable from 
those of perfect fluids. Let us reemphasize that we do not propose elastic matter 
as a physically realistic matter model for the universe; however, the results bear 
relevance on our general understanding of the physics of cosmological models. 
The isotropization rates of elastic cosmological models differ from those of per- 
fect fluid models even in the case where the material properties of the elastic 
matter deviate from those of perfect fluid matter by an arbitrarily small amount. 
Whether the difference in isotropization rates found here carries over to more 
general cosmological models with more general matter sources remains to be 
seen. 

Acknowledgments: We would like to thank an anonymous referee for valuable 
comments. S. C. is supported by Ministerio Ciencia e Innovation, Spain (Project 
MTM2008-05271). 

Appendix: Elasticity theory 

This appendix is devoted to a presentation of some basic concepts of the general 
relativistic theory of elasticity. In particular, we specify in detail the assump- 
tions on the elastic matter model used in this paper, which lead to the energy 
density © and the principal pressures ([7]) . A more detailed exposition of rela- 
tivistic elasticity can be found in the references listed at the end of the paper. 

The reference state of an elastic body is defined to be the state of the body 
in the absence of strain and external forces. (The reference state has of course 
to be understood in a Platonic sense, since the conditions of vanishing strain 
and external forces cannot be realized in the the real world.) It is assumed that, 
in the reference state, the particles of the body form, in the continuum limit, a 
smooth three dimensional manifold N, called material space. The manifold N 
must be equipped with different tensor fields in order to describe the physical 
properties of the body in the reference state. The very least one has to require 
is the ability to "count" the particles, and therefore to admit that the material 
space is equipped with a volume form (particle density). However, a consistent 
theory for the dynamics of clastic bodies in general relativity is presently avail- 
able only under the stronger requirement that on the manifold N be defined 
a Riemannian metric 7 (which reflects the ability to measure the "distance" 
between the particles). We refer to [TO] for a discussion of the physical inter- 
pretation of this condition and on the restrictions that it imposes on the class 
of elastic materials covered by the theory. 

The state of the body in a four dimensional space-time (M, g) is determined 
by a configuration map tp : M —* N with the property that for all q € N the set 
tp~ 1 (g) is a timelike curve (the world-line of the "particle" q). This definition 
implies that the kernel of the deformation gradient Tip : TM — > TN is generated 
by a (future-directed unit) timelike vector field it, which is interpreted as the 
matter four-velocity; by construction, ip _1 (<?) is an integral curve. 

In the following, let x 11 denote a system of local coordinates in the space-time 
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and X A a system of local coordinates in the material space; Greek indexes run 
from to 3, capital Latin indexes from 1 to 3 and refer to tensors on the material 
space. Let ip A (x 11 ) be the triple of functions representing the configuration 
map in these coordinates and d^^ix^) the matrix- valued function representing 
the deformation gradient. Moreover we denote by jab the components of the 
material metric. We define two metrics on the orthogonal complement (u) of 
u in TM. The Riemannian metric induced by g we denote by g: 

Q^iv g ' ^iv U^XLi/ . 

The pull-back of the material metric by the map ip, i.e., ip*^), is called the 
relativistic strain tensor h: 

K„ = d^ A d v ip B 1AB • (38) 
The relativistic strain tensor satisfies two fundamental properties: 
(i) It is orthogonal to the matter four-velocity, i.e., 

h^u" = ; (39a) 



(ii) It is constant along the matter flow, i.e., 

Cuh MU = . (39b) 



By the property (i), defines a Riemannian metric on (u)' L . Hence h% = 
g^hau has three positive eigenvalues hi, h 2 , /13. 

The material is unstrained at x iff g^ix) — h^vix). The scalar quantity 

n = \J det g h — \fh\h2h3 

is the particle density of the material. This interpretation is justified by virtue 
of the continuity equation 

V p (m") = . 

A specific choice of elastic material is made by postulating a constitutive 
equation, i.e., the functional dependence of the (rest frame) energy density p of 
the material on the configuration map, the deformation gradient and the space- 
time metric. An important class of materials is the one for which this functional 
dependence enters only through the principal invariants of the strain tensor. In 
this case we have 

p = p(gi,?2,?3) , (40) 

where 

qi =trh, q 2 = tr (h 2 ) , q 3 = tr (h 3 ) ; 

since n 2 = (qf — 3qiq2 + 2q 3 )/6, one of the invariants qi can be replaced by 
the particle density n. The materials described by (|40]) generalize the class 
of isotropic, homogeneous, hyperelastic materials from the classical theory of 
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elasticity, see [TO] . As explained in [2] , restriction to these materials guarantees 
that the resulting elasticity theory is generally covariant and therefore seems 
quite a natural assumption in the context of general relativity. 

In this paper we make use of a constitutive equation introduced in [TO] . Let 
the shear scalar to be defined as 



= ^K 2 (<?i-93) -24] , 



(41a) 



or, in terms of the eigenvalues hi, hi, ha 



1 

12 




(41b) 



Evidently, s 2 is non-negative, and s 2 = (no shear) iff oc g^ (or equiva- 
lently, hi = h-x = /13). Following [TO] we shall consider a constitutive equation 
of the quasi-Hookean form, that is 



P = P( n ) + fi(n)s 2 , 



(42) 



where p(n) is the unsheared energy density and fi{ri) the modulus of rigidity 
(or shear modulus). The stress-energy tensor associated with these materials is 
obtained as the variation with respect to the space-time metric of the matter 
action Sm — — J \fW\P- The result is given in [TOJ Sec. 6] and reads 



where = p g^ v 



(43a) 



6 n 2 



(tr(/t 3 ) - (tr/i) 3 ) 9liu + {tihfh^ - (h 3 )^ 



Here p is the isotropic (component of the) pressure, which is given by 



p = p{n) + i>{n)s 2 



2 d ( p\ ( dfi 

where p = n — 
an 



(43b) 



u=[n--H). (44) 



The principal pressures Pi (which are the [non-zero] eigenvalues of T^,) are thus 
of the form Pi = p + 5pi, for an unstrained configuration, pi — p, i — 1, 2,3. 
For jl — (or s 2 = 0), the elastic material reduces to a perfect fluid with 
stress-energy tensor T^ v — pu^u^ + pg^, energy density p — p and pressure 
p = p. 

It remains to specify the functions p and p, in the constitutive equation (|42p . 
Following [TO] , we postulate a linear equation of state between the modulus of 
rigidity p. and the unsheared pressure p, 



p = vp . 



(45) 
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where v is a dimensionless constant that we call elastic constant and that is al- 
lowed to vanish (in which case the elastic body becomes a perfect fluid) . Finally 
we postulate a linear equation of state between the unsheared pressure p and 
the unsheared energy density p, 

p = wp , (46) 

where the constant w is allowed to take values in the interval [—1,1]. We also 
assume that the product vw is non-negative to guarantee the non-negativity of 
the energy density p (see Eq. (|4T)) below) . We remark that our choice of equation 
of state (|4"6"|) is different from the one in [TU] ; in ref . [TU] , the unsheared pressure 
and the unsheared energy density are assumed to satisfy a polytropic equation 
of state. The equation of state p = p(p) selects the perfect fluid model that 
arises as a special case of the elastic matter model when the modulus of rigidity 
vanishes. In our case, it is a perfect fluid that obeys a linear equation of state, 
which is the perfect fluid model most widely used in cosmology. In [10] , where 
the applications concern the stellar dynamics case, the elastic body becomes a 
polytropic perfect fluid star in the absence of rigidity. 
Substituting (@SJ| and IpEg] ) in (|i3] l we find 

p = p Q n w+1 , p = p vwn w+1 (\w\ < 1, vw > 0) 

for some constant po > 0. Accordingly, 

p = p n w+1 (l + vw s 2 ) , p = wp. (47) 

Since for an unstrained material p = p and pi—p — p hold, i = 1,2, 3, the bound 
|w| < 1 ensures that the dominant energy condition \pi\ < p is satisfied for an 
unstrained configuration. Furthermore, vw > guarantees that the energy 
density is positive for all values of the shear scalar s 2 and has a minimum at 
zero shear. When v = 0, the modulus of rigidity fi vanishes and the elastic 
matter reduces to a perfect fluid with a linear equation of state p — wp; the 
condition \w\ < 1 ensures that the dominant energy condition \p\ < p is satisfied 
for this perfect fluid. When w — (so that p = 0), the choice of v is irrelevant, 
since vw — 0; this is clear because shear cannot occur for dust. Henceforth, 
unless stated otherwise, by elastic matter we will always mean matter with 
constitutive equation ([6|), where w e [—1,1] and vw > 0. 

Consider now a homogeneous space-time (M,g) of Bianchi type I, i.e., 

g^dx^dx" = -dt 2 + g l:j (t)dx l dx 3 , (48) 

where gij{t), i,j — 1,2,3, is a family of Riemannian metrics that is induced 
on the spatially homogeneous hypersurfaces t = const. The spatial coordinates 
are chosen so that the Killing vector fields of the space-time coincide with the 
operators d x i . Now we recall that a matter model in a spatially homogeneous 
space-time is said to be non-tilted if the matter four-velocity u is orthogonal to 
the hypersurfaces of spatial homogeneity. For the metric (|4"8|) this means that 
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u = dt- In this paper we assume that the elastic material in the Bianchi type I 
space-time is non-tiltcd. By (l39l) this implies 

d t V = , h „ = , (49) 

in particular the relativistic strain tensor is time independent. Next we recall 
(see, e.g., [TB]) that compatibility with the Einstein equations entails that a 
matter model with stress-energy tensor in a Bianchi type I space-time must 
satisfy d x iT^ v = 0, i.e. the stress- energy tensor is independent of the spatial 
variables. A matter model is said to inherit the Bianchi type I symmetry if the 
relation d^T^ = implies that the matter fields, which enter into the definition 
of T^ v , are also independent of the spatial variables. For instance, in the case of 
perfect fluids, the relation d x iT^ v = implies the relations d x tp — d x ip = on 
the energy density and the pressure. In this paper we restrict to elastic matter 
that inherits the Bianchi type I symmetry, in particular we assume that 

9 X »V=0. (50) 

The equations (|49|) and (|50f imply that the component of the relativistic strain 
tensor are constant: 

hoa = hok = 0, hj = const . (51) 

These constants can be fixed arbitrarily and correspond to the "initial data" 
for the matter field. Note that the relativistic strain tensor has been fixed by 
our geometric assumptions and that the concept of material space and material 
metric become redundant. (This is of course a consequence of the high degree of 
symmetry of the space-time). However at this point it is instructive to derive a 
material metric and configuration maps that give rise to a relativistic strain ten- 
sor of the form (f5"Tj) . We may consider a flat material metric and fix coordinates 
on the material space such that jab = Sab', then we restrict to homogeneous 
deformations (see |12j ) 

d t <ip A = , d x iip A = F t A = const =*> 4j a = F i A x l + c A , 

where c A is a constant. By ([38]) we obtain the relativistic strain tensor (fSTj) , 
where % = SabF^F^. Substituting (|5l"|) in ([43]) we obtain 

Too = P, T k = jk = 0, fij = T %j (52) 

where T; j is given in terms of hij via (|43b|) . 

The results presented in this paper have been obtained for diagonal solutions 
of the Einstein-elastic matter equations, i.e., solutions for which the space-time 
metric and the stress-energy tensor are diagonal. These solutions correspond 
to initial data, say at time t = 0, of the following type. Let gij(0) be the ini- 
tial spatial metric and kij (0) the second fundamental form of the hypersurface 
t = 0. Without loss of generality we can assume that 3^(0) and k l j(0) are diag- 
onal (by choosing coordinates adapted to an orthogonal basis of eigenvectors of 
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£^(0)). Furthermore we impose the condition that hij is diagonal; in particu- 
lar, by rescaling the spatial coordinates, we can assume hij = <5y ;. The Einstein 
equations, in units c = 1 = 87rG, decompose into the momentum constraint 
jk = 0, which is automatically satisfied by (|52p , the Hamiltonian constraint 

(trfc) 2 - k l 3 k\ - 2p = , (53a) 

and the evolution equations 

d t9ij = ~2k l3 d t k\ 3 = (trfc)fc l J - T) + lg l j {T\ - p) . (53b) 

Since the off-diagonal elements of the tensor T* • form an homogeneous polyno- 
mial in h 1 ^ = g lk hjk, i ^ j, it follows from the evolution equations (|53b[) that 
[gij, k 1 ^, h l j) (and therefore T 1 -) remain diagonal for all times. We refer to this 
special class of solutions of the equations (f53|) as diagonal models. 

From h l = g tk hkj = diag(<7 n , g 22 , g 33 ) = diag(/ii, ti2, /13) we conclude that 
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(54) 



cf. (|41b[) . which can be inserted into (|4"7} . i.e., 

P = P o(g 11 g 22 g 33 ) {w+1)/2 (i + vw.s 2 ), (55a) 

to yield p as a function of g 11 , g 22 , g 33 . Moreover, from (I43b[) we find 

1 /^ii ff 33 g ll g 22\ 

T i=P+gA^- ^11 + ^2 -gil)> (55b) 

+ (55c) 
2_ 6 U 11 g 22 g 33 g 22 ) 1 



T 3=P+gM^-^33+^TT-^33j > ( 55d ) 

where p = u>p and /i = pov w{g 11 g 22 g 33 )^ w+1 ^ 2 and are thus functions of g 11 , 

g 22 ,g 33 - 

Finally, in the local rotational symmetry case, where we assume that the 
plane of rotational symmetry is the X2-X3 plane, the metric takes the form ((T|). 
Substituting g 11 — 1/A and g 22 = g 33 = l/B in Eq. (|55a|) we obtain the 
expression © for the energy density. Furthermore, defining pa = T \, ps = T 2 2 , 
the rescaled principal pressures wa — Pa/ P, hb = Pb/p are given by (|7|). 
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